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Recent experiments [I.V. Grigorieva et al, Phys. Rev. Lett. 96, 077005 (2006)] on visualization of 
vortices using the Bitter decoration technique revealed vortex shells in mesoscopic superconducting 
Nb disks containing up to L = 40 vortices. Some of the found configurations did not agree with those 
predicted theoretically. We show here that this discrepancy can be traced back to the larger disks 
with radii ii ~ 1 to 2.5pim, i.e., R ~ 50 — 100^(0) used in the experiment, while in previous theoretical 
studies vortex states with vorticity L < 40 were analyzed for smaller disks with i? 5 — 20^(0). The 
present analysis is done for thin disks (mesoscopic regime) and for thick (macroscopic) disks where 
the London screening is taken into account. We found that the radius of the superconducting disk 
has a pronounced influence on the vortex configuration in contrast to, e.g., the case of parabolic 
confined charged particles. The missing vortex configurations and the region of their stability are 
found, which are in agreement with those observed in the experiment. 
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I. INTRODUCTION 

A mesoscopic superconducting disk is the most simple 
system to study confined vortex matter where effects of 
the sample boundary plays a crucial role. At the same 
time, it is a unique system because just by using disks 
of different radii, or by changing the external parame- 
ters, i.e., the applied magnetic field or temperature, one 
can cover — within the same geometry — a wide range 
of very different regimes of vortex matter in mesoscopic 
superconductors. 

Early studies of vortex matter in mesoscopic disks were 
focused on a limiting case of thin disks or disks with 
small radii in which vortices arrange themselves in rings 
P, 0, H, 13, H , in contrast to infinitely extended super- 
conductors where the triangular Abrikosov vortex lat- 
tice is energetically favorable 0, H, Several stud- 
ies were devoted to the questions: i) how the vortices 
are distributed in disks, ii) which vortex configuration 
is energetically most favorable, and iii) how the tran- 
sition between different vortex states occurs. Lozovik 
and Rakoch Q analyzed the formation and melting of 
two-dimesional microclusters of particles with logarith- 
mic repulsive interaction, confined by a parabolic poten- 
tial. The model was applied, in particular, to describe 
the behaviour of vortices in small thin (i.e., with a thick- 
ness smaller than the coherence length £) g rains of type 
II superconductor. Buzdin and Brison [lO[ studied vor- 
tex structures in superconducting disks using the image 
method, where vortices are considered as point-like "par- 
ticles", i.e., within the London approximation. Palacios 
[m calculated the vortex configurations in superconduct- 
ing mesoscopic disks with radius equal to R — 8.0^, 
where two vortex shells can become stable. The de- 
magnetization effects were included approximately by in- 
troducing an effective magnetic field. Geim et al. [l^ 
studied experimentally and theoretically the magnetiza- 



tion of different vortex configurations in superconducting 
disks. They found clear signatures of first- and second- 
order transitions between states of the same vorticity. 
Schweigert and Peeters Q analyzed the transitions be- 
tween different vortex states of thin mesoscopic super- 
conducting disks and rings using the nonlinear Ginzburg- 
Landau (GL) functional. They showed that such transi- 
tions correspond to saddle points in the free energy: in 
small disks and rings — a saddle point between two gi- 
ant vortex (GV) states, and in larger systems — a saddle 
point between a multivortex state (MV) and a GV and 
between two MVs. The shape and the height of the nucle- 
ation barrier was investigated for different disk and ring 
configurations. Milosevic, Yampolskii, and Peeters 
studied vortex distributions in mesoscopic superconduct- 
ing disks in an inhomogeneous applied magnetic field, 
created by a magnetic dot placed on top of the disk. It 
was shown 14], that such an inhomogeneous field can 
lead to the appearance of Wigner molecules of vortices 
and antivortices in the disk. 

In the work of Baelus et al. 15] the distribution of 
vortices over different vortex shells in mesoscopic super- 
conducting disks was investigated in the framework of 
the nonlinear GL theory and the London theory. They 
found vortex shells and combination of GV and vortex 
shells for different vorticities L. 

Very recently, the first direct observation of rings of 
vortices in mesoscopic Nb disks was done by Grigorieva 
et al. [l^ using the Bitter decoration technique. The 
formation of concentric shells of vortices was studied for 
a broad range of vorticities L. From images obtained 
for disks of different sizes in a range of magnetic fields, 
the authors of Ref. [l^ traced the evolution of vortex 
states and identified stable and metastable configurations 
of interacting vortices. Furthermore, the analysis of shell 
filling with increasing L allowed them to identify magic 
number configurations corresponding to the appearance 



2 



of consecutive new shells. Thus, it was found that for 
vorticities up to L = 5 all the vortices are arranged in 
a single shell. Second shell appears at i = 6 in the 
form of one vortex in the center and five in the second 
shell [state (1,5)], and the configurations with one vortex 
in the center remain stable until i = 8 is reached, i.e., 

(1.7) . The inner shell starts to grow at L = 9, with the 
next two states having 2 vortices in the center, (2,7) and 

(2.8) , and so on. From the results of the experiment [ll| 
it is clear that, despite the presence of pinning, vortices 
generally form circular configurations as expected for a 
disk geometry, i.e., the effect of the confinement dom- 
inates over the pinning. Similar shell structures were 
found earlier in different systems as vortices in super- 
fluid He [13, E, H H 111, charged particles confined 
by a parabolic potential [22|, dusty plasma [11], and col- 
loidal particles confined to a disk [24]. Note that the 
behavior of these systems is similar to that of vortices in 
thin mesoscopic disks, thus our approach of Sec. II can 
be used for better understanding of the behavior of vari- 
ous systems of particles confined by a parabolic potential 
and charachterized by a logarithmic interparticle interac- 
tion (e.g^vortices in a rotating vessel with superfluid He 
[itI. [isl. TiqI. [20|. [2]| ) . In contrast, our results presented in 
Sec. HI are specific for vortices in thick large mesoscopic 
superconducting disks, where the London screening is im- 
portant, and the intervortex interaction force is described 
by the modified Bessel function. 

The filling of vortex shells was experimentally ana- 
lyzed [l^ for vorticities up to L = 40. Many configura- 
tions found experimentally agree with earlier numerical 
simulations for small L which were done for mesoscopic 
disks with radii as small as R ~ 6 — 8^(0), although 
the disks used in the experiments [l6| were much larger, 
i? w 50 — 100^(0). At the same time, some theoretically 
predicted configurations were not found in the experi- 
ment, such as states (1,8) for L = 9 and (1,9) for L — 10. 
The difference between vortex states in small and large 
disks becomes even more striking for larger vorticities L. 
In small disks with radii of a few ^, the formation of GVs 
is possible if the vorticity L is large enough, e.g., in disks 
with R = 6^, a GV with L = 2 appears in the center 
for total vorticity L = 14, but for vorticities L > 14, all 



the vortices form a GV [15[. Obviously, this boundary- 



induced formation of GVs is possible only in the case of 
small disks: in large disks vortices instead form the usual 
Abrikosov lattice which is distorted near boundaries. 

The aim of the present paper is to theoretically analyze 
vortex states, using Molecular-Dynamics (MD) simula- 
tions, in rather large mesoscopic superconducting disks 
and thus to study the crossover between mesoscopic and 
macroscopic disks, i.e., the regime corresponding to the 
Nb disks used in the recent experiments of Ref. Ja] and 
to look for the missing vortex configurations in the earlier 
simulations. We analyzed the region of stability of those 
configurations and performed a systematic study of all 
possible vortex configurations. We found that the radius 
of the disk has an influence on the vortex shell struc- 



ture, in contrast to the case of charged particles confined 
by a parabolic potential [2^ . This analysis was done 
for thin mesoscopic disks and for thick disks where the 
London screening becomes pronounced. We also perform 
calculations of vortex configurations using the GL equa- 
tions, and we compare these results to the ones obtained 
within the MD simulations. The calculated vortex con- 
figurations agree with those observed in the experiment 

The paper is organized as follows. In Sec. II, we dis- 
cuss thin mesoscopic disks. The model and the sim- 
ulation method is described in Sec. II. A. In Sec. II. B, 
we discuss different vortex configurations and the forma- 
tion of vortex shells. We analyze the ground states and 
metastable states for different vorticities in Sec. II. C, us- 
ing a statistical study, similar to the one employed in the 
experiments of Ref. [Ifl , by starting with many random 
vortex configurations and comparing the energies of dif- 
ferent vortex configurations. Based on that analysis we 
reconstructed the "radius i?— magnetic field H" phase 
boundary (Sec. II. D). In Sec. HI, we study the crossover 
from thin mesoscopic to thick macroscopic disks, and we 
analyze the impact of the London screening on the vor- 
tex patterns. In Sec. IV, we calculate the crucial vortex 
configurations in disks using the GL equations, and we 
compare them to the results obtained within the MD sim- 
ulations. A summary of the results obtained in this work 
is given in Sec. V. 



II. MESOSCOPIC DISKS 

A. Theory and simulation 

In this Section, we consider a thin disk with thickness 
d and radius R such that Ae// ^ R ^ £, ^ d, placed 
in a perpendicular external magnetic field Hq. Here 
Ae// = X^/d is the effective London penetration depth 
for a thin film, A is the bulk London penetration depth, 
and ^ is the coherence length. We follow here the theoret- 
ical approach developed in [l^, [l^, |2^ for thin disks and 
we use the original dimensionless variables used in those 
works. Thus, following [H, [1^ the lengths are measured 
in units of the coherence length ^, the magnetic field in 
units of Hc2 = ch/2e$,^ = kV^Hc, and the energy density 
in units of H'^/8tt. The number of vortices, or vorticity, 
will be denoted by L. In a thin disk in which demag- 
netization effects can be neglected the free energy in the 
London limit can be expressed as [l^, [l^, [2^ 
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where the potential energy of vortex confinement consists 
of two terms: 



s self . shield 
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i.e., the interaction energy between the ith vortex and 
the radial boundary of the superconductor, 



self 



(3) 



and the interaction energy between the ith vortex and 
the shielding currents. 



^shield 



In Eq. HD, 
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is the repulsive interaction energy between vortices i and 
j. Here = pi/R is the distance to the vortex normal- 
ized to the disk radius. The divergence arising when i = j 
is removed in Eq. (O using a cutoff procedure (see, e.g., 
SEIli) which assumes the replacement oi \ pi — pj \ 
by a (or by in not normalized units) for i ^ j. Finally, 
gcore ^ (2/RfL ln{R/a) and e^*^''^ = R^HH'^ are the 
energies associated with the vortex cores and the exter- 
nal magnetic field, respectively. Notice that the energy 
of the vortex cores e'=°'"'^ becomes finite due to the cutoff 
procedure and is strongly dependent on the cutoff value 
[l^. Here we use for the vortex size a — -\/2^; this 
choice, as shown in Ref. [1^ , makes the London and the 
Ginzburg-Landau free energies to agree with each other. 

From the expression of the free energy given by 
Eqs. inn]), we obtain the force acting on each vortex, 
— VfeG(pi,Pj), where — V/t is the gradient with respect 
to the coordinate p^. This yields a force per unit length. 
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in units of i7^^/87r, where the summation over k runs 
from 1 to L, except for k — i. The first term describes 
the vortex interaction with the current induced by the 
external field and with the interface, 
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The second term is the vortex- vortex interaction 
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The above equations allow us to treat the vortices as 
point-like particles and the forces resemble those of a 
two-dimensional system of charged particles with 1/r re- 
pulsive interaction confined by some (usually, parabolic) 
potential [2^ . However, the inter- vortex interaction in 
our system is different from 1/r, and the confined po- 
tential differs from parabolic and depends on the applied 
magnetic field. 



To investigate different vortex configurations, we per- 
form MD simulations of interacting vortices in a disk (see, 
e.g., Ref. [IBl), starting from randomly distributed vor- 
tex positions. The final configurations were found after 
typically 10^ MD steps. 

In order to find the ground state (or a state with the 
energy very close to it) we perform many (typically, one 
hundred) runs of simulations for the same number L of 
vortices starting each time from a different random ini- 
tial distribution of vortices. As a result, we obtain a 
set of final configurations which we analyze statistically, 
i.e., we count probabilities to find the different configura- 
tions with the same vorticity L, e.g., configurations (1,8) 
and (2,7) for L = 9. We can expect that the configu- 
ration which appears with the highest probability is the 
ground state of the system, i.e., the vortex state with the 
lowest energy. (However, in some cases, i.e., for partic- 
ular vortex configurations the highest probability state 
turns out to be not always the ground state configura- 
tion. One of these special cases will be addressed below.) 
This approach corresponds to the analysis done in the 
experiment [16]. 

The MD simulation was performed by using the 
Bardeen-Stephen equation of motion 



dpi^ 
dt 



(9) 



where i denotes the ith vortex, 7] is the viscosity coeffi- 
cient rj ^ ^oHc2/ PnC^ , with p„ being the normal-state 
resistivity; $o = hc/2e is the magnetic flux quantum. 
The time integration was accomplished by using the Eu- 
ler method. 



B. Vortex configurations for different L: 
of vortex shells 



Formation 



To study the formation of vortex shells in mesoscopic 
supeconducting disks, here we analyze the evolution of 
vortex configurations with increasing number of vortices, 
L, in a disk with radius R = 50^. The results of our cal- 
culations for L = 1 to 10 are presented in Fig. 1. When 
the vorticity L of the sample increases, the vortex con- 
figurations evolve with increasing applied magnetic field 
as follows: starting from a Meissner state without vor- 
tex, then one appears in the center (Fig. 1(a)), which we 
denote as (1), then two symmetrically distributed vor- 
tices, (2) (Fig. 1(b)). Further increasing the magnetic 
field results in the formation of triangular, (3) (Fig. 1(c)) 
and square like (4) (Fig. 1(d)) vortex patterns in the sam- 
ple, and in a five- fold symmetric pattern, (5), shown in 
Fig. 1(e). When the vorticity L increases from 5 to 6, 
a vortex appears in the center of the disk, thus starting 
to form a second shell of vortices in the disk (Fig. 1(f)). 
We denote the corresponding two-shell vortex state con- 
taining 1 vortex in the first shell and 5 vortices in the 
second shell as (1,5). Two-shell configurations with one 
vortex in the center and newly generated vortices added 
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FIG. 1: The evolution of vortex configurations for the states 
with vorticity increasing from L = 1 to 10, in mesoscopic 
superconducting disks with radius R = 50^. Vortices form 
one-shell configurations for L = 1 to 5 (a)-(e). The formation 
of second shell starts for L = 6 (f). The inner shell contains 
1 vortex at the center for L = 5 to 8, while new vortices fill 
the outer shell (f-h). For L — 9, second vortex appears in the 
inner shell (i), and then second shell again starts to grow for 
L = 10 (j). 



to the outer shell remain for the states (1,6) with L = 7 
(Fig. 1(g)), and (1,7) with L = 8 (Fig. 1(d)). The num- 
ber of vortices in the inner shell begins to grow at L = 9 
thus forming subsequently configurations (2,7)(Fig. l(i)), 
and (2,8) with L = 10 (Fig. l(j)). 

Note that in earlier theoretical works on vortices in 
mesoscopic superconducting_^disks, a configuration (1,8) 
for L = 9 was predicted |J, [l^ as a ground state in 
smaller disks, which however was not observed in the ex- 
periment 'l^l as a stable state (the special case L = 9 
will be discussed in Sec. Ill for the model of a thick disk, 
i.e., d ^ X, relevant to the experiment p^ . where the 
London screening in large disks, i.e., i? > A, is taken into 
account). Our calculations show that the multivortex 
states with two vortices in the center and the other vor- 
tices on the outer shell can exist till L ~ 14. For L > 14, 
the inner shell starts growing again till L = 16, which 
means that a newly nucleated vortex will be generated 
in the inner shell, while the number of vortices in the 
outer shell stays the same. We found that those config- 
urations are (3,11) for L = 14, (4,11) for L = 15, and 
(5,11) for L — 16. At L — 17, a vortex appears in the 
center, thus giving rise to the formation of a third shell 
with one vortex in the center. With increasing number 
of vortices in the disk, the next three vortices are added 
to the outermost shell, after which all three shells grow 
intermittently till L = 32. The fourth shell appears at 
L = 33 in the form of one vortex in the center. The bor- 
derline vortex configurations illustrating the formation 
of new shells are presented in Fig. 2. We summarize the 
vortex configurations found for vorticities L = 1 to 33 in 
Table 1. 

It is appropriate to mention that in our numerical 
calculations using the vortex-vortex interaction force, 
Eq. ([5]), the obtained vortex patterns for some L are dif- 
ferent from those found in Ref. [3] for particles with log- 
arithmic interaction, confined by a parabolic potential, 
even in case when the interaction with images was taken 
into account. Although for many vorticities L both ap- 
proaches result in the same "robust" configurations (i.e., 
less sensitive to details of interparticle interactions) , there 
is essential difference in configurations for some other vor- 
ticities. The lowest vorticity for which our results devi- 
ate from those obtained in Ref. [Sj] is L = 9: this spe- 
cial case is described in detail below. Note, for instance, 
differences in three-shell configurations: (1,5,12) (our ap- 
proach, see Table 1) [and (1,6,11) Q] for L = 18, (1,7,12) 
[(1,6,13)] for L = 20, (1,8,13) [(1,7,14)] for L = 22, 
(2,8,13) [(1,8,14)] for L = 23, (3,10,15) [(4,9,15)] for 
L = 28, (5,10,16) [(4,10,17)] for L = 31, and (5,11,16) 
[(4,11,17)] for L = 32. Moreover, the filling of next 
(fourth) shell starts, according to our calculations, for 
i = 33 (i.e., (1,5,11,16)) while in Ref. [3] this transition 
occurs for L = 34. As noticed in Ref. [3], the interaction 
with images leads to the stabilization of configurations 
with larger number of vortices on inner shells. While in 
Q this tendency was revealed only in rather large clus- 
ters (i.e., the appearance of different configurations for 
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FIG. 2: The borderline vortex configurations for: one shell, 
starting from state (1) for L = 1 to state (5) for L — 5 (a), 
two shells, states (1,5) for L — 6 (b) to (5,11) for L = 16 
(c), three shells, states (1,5,11) for L = 17 (d) to (5,11,16) 
for L — 32 (e), four shells, starting from state (1,5,11,16) for 
L — 33 (f), for disks with radius R = 50^. 

L — 45: i) (2,8,14,21), without interaction with images, 
and ii) (3,8,14,20), if the interaction with images is taken 
into account), we found that the stabilization of vortex 
shell structures with a larger number of vortices on inner 
shells occurs for vorticities L > 20, due to the interaction 
with images and with the shielding current induced by 
increasing magnetic field at the boundaries of the disk 

da. 

Note that, although here we employ the model of a thin 
mesoscopic disk, the results of our calculations for the fill- 
ing of vortex shells perfectly match those discovered in 
the experiment ^1^ (where the disks were rather thick) 
for vorticities L = 1 to 40. The vortex configurations 
calculated in this section using many runs of simulations 
with random initial distributions, as will be shown be- 
low, are not always the ground-state configurations. The 
obtained results imply that the size of the disk influences 
the vortex configurations in superconducting disks. This 
will be clearly demonstrated in Sec. Ill where we consider 
thick disks and we show that the London screening has a 
pronounced impact on the vortex configurations in disks 



of different radii. 



C. The ground state and metastable states 

As described in Sec. II. A, in order to find the ground 
state of the system (or a state with energy very close 
to it), we performed many (usually one hundred) simu- 
lations for the same number of vortices. In most cases, 
always one configuration dominated over the other pos- 
sible configurations for a certain vorticity L, which was 
identified as the "ground state" . However, for some vor- 
tex configurations competing states appeared with com- 
parable probabilities. 

Let us now consider those special cases. For in- 
stance, it follows from our calculations that two con- 
figurations, (1,8) and (2,7), are possible for the same 
vorticity L ^ 9. They are shown in Figs. 3(c), (d) (an- 
other example of competing states with the same vortic- 
ity, e.g., L — 17, are the three-shell configuration (1,5,11) 
and the two-shell configuration (5,12) which are shown 
in Figs. 3(a), (b), correspondingly). We found that the 
configuration (2,7) is the ground state for L = 9 in a disk 
with radius i? < 11^^ (see the R — H phase diagram in 
Fig. 6) in a certain range of magnetic fields. In very large 
disks with R > 50^ the vortex state (2,7), although being 
a metastable state, is the highest probability state. Note 
that configuration (2,7) was also found as a ground state 
for the system of charged particles [22| . At the same time, 
as it was shown in Ref. |15|] using the GL theory, in small 
superconducting disks (e.g., for radius R = 6^), this con- 
figuration occurred to be a metastable state, while state 
(1,8) was found as the ground energy state. This clearly 
shows that the ground state configuration for a certain L 
depends on the radius of the disk. 



1. Statistical study of different vortex states 

In the experiments [T6| . vortex configurations were 
monitored in large arrays of similar mesoscopic disks 
(dots). This allowed them to study the statistics of the 
appearance of different vortex configurations in practi- 
cally the same disk. The results show that, e.g., in a disk 
with radius R = 1.5/iTO and magnetic field Hq = 60Oe, 
configuration (2,8) for L = 10 appeared more frequently. 
Other configurations for the same total vorticity L = 10, 
e.g., configuration (3,7) appeared only in few cases. In- 
terestingly, not only various configurations with the same 
total vorticity L = 10 appeared, but also vortex states 
with L = 9 (2,7) and — less frequently — two modifica- 
tions of the state (1,8): with a ring-like outer shell as well 
as a square-lattice-like vortex pattern. This statistical 
study provides indirect information about the ground- 
state and metastable states. 

We performed a similar investigation of the statistics of 
the appearance of different vortex states for ideal disks, 
i.e., in the absence of pinning. One hundred randomly 
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(2,8,13) 


(2,8,14) 
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(3,9,14) 
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TABLE I; Formation of vortex shells in mesoscopic superconducting disks: vortex configurations for different vorticities L. 




FIG. 3: Possible vortex configurations for the total vorticity 
L = 17 in a disk with radius R = 50^: (1,5,11) (a) and (5,12) 
(b); and for the total vorticity L = 9: (2,7) (c) and (1,8) (d). 

distributed initial states were generated for our statis- 
tical study for each set of parameters. We studied the 
dependence of the appearance of different vortex config- 
urations on the applied magnetic field. For instance, for a 
disk with radius R — 42^ and the magnetic field varying 
from H — 0.011i/c2 to 0.015i/c2, we counted how often 
the different configurations (e.g., (1,8) and (2,7)) for a 
total vorticity L = 9 appeared. 

The results of such calculations are shown in Fig. 4. 
At low magnetic field, H = 0.011i?c2, the disk cannot 
accommodate 9 vortices, so the number of configurations 
(1,8) and (2,7) is zero, and in most cases we obtain the 
configurations (1,6) or (1,7) for L = 7 and L = 8, re- 
spectively. As the magnetic field increases, the probabil- 
ity to find the configurations (1,8) and (2,7) increases, 
and at the same time the probability to find the config- 
urations (1,6) and (1,7) decreases. Our statistical result 
shows clearly that in the range of magnetic field shown in 
Fig. 4 the probability to find configuration (2,7) is always 
higher than to find configuration (1,8). 

Similar analysis was performed for a disk with radius 
R = 47^ where we looked for configurations with total 
vorticity L = 17 in the range of the applied magnetic 
fields from H = 0.014i?c2 to 0.023Hc2- As shown in 
Fig. 5, configurations (4, 10) for L = 14, (5, 10) for L = 
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FIG. 4: The results of the statistical study of different vortex 
configurations, i.e., the probability to find a given vortex state 
as a function of the appfied magnetic field, for a disk with 
radius R = 42^ and varying magnetic field. 
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FIG. 5: The results of the statistical study of different vortex 
configurations, for a disk with radius R = 47^ and varying 
magnetic field. 

15, (5, 11) for L = 16, and (1, 5, 11) and (5, 12) for L 17 
dominate with increasing magnetic field. Note that two 
configurations, (1,5,11) and (5,12), appeared in the same 
magnetic field range for L = 17, and (1,5,11) is always the 
dominant configuration, i.e., the formation of the third 
shell starts for vorticity L = 17 (cp. Ref . Kdi) . The results 
of the statistical study of the configurations (2,7) and 
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FIG. 6: The results of the statistical study of different vor- 
tex configurations, for a disk with radius i? = 9^ and varying 
magnetic field. The total probabilities of vortex states with 
different vorticities L = 7, 8, 9, and 10, are plotted. For 
vorticity L = 9, the probabilities of the two possible configu- 
rations (2,7) and (1,8) are also shown. 

(1,8) for small disks (i? = 9^) are shown in Fig. 6. 
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FIG. 7: The R — H phase diagram for states with total vor- 
ticity L = 9, for small disks with radii changing in the range: 
R = 7^ to R = 12^. Points A and B correspond to the config- 
urations (a) and (b) shown in Fig. 10, respectively. The inset 
shows the confinement potentials (see Eq. ([2])) corresponding 
to the points C, D, and E in the phase diagram. 



2. The R — H phase diagram 

To find the region of the existence and stability of vor- 
tex states with different vorticity L, we performed a di- 
rect calculation of their energies using Eq. ^ . 

As an example, we considered the configurations (1,8) 
and (2,7) in disks with radius changing in a very wide 
range from i? = 4^ to i? ~ 100^. We change the radius 
of the disk, and at the same time keep the flux passing 
through the disk $ — SHq the same, in order to keep the 
same vorticity L in the disk. Here <I> is the flux passing 
through the specimen, Hq is the applied magnetic fleld, 
S — ttRq is the surface area of the specimen. 

The phase diagram "radius of the disk R - applied 
magnetic fleld H" is shown in Fig. 7 for i? = 7^ to 12^. 
According to our calculations, for small radii R < 7^, the 
energy of the conflguration (2,7) appears to be always 
lower than that of conflguration (1,8). The total energy 
for both conflgurations, (2,7) and (1,8), decreases with 
increasing magnetic fleld, and the energy of the state 
(2,7) is slightly larger than that of the state (1,8). 

For disks with radius between 7^ and 12^, the conflg- 
uration (1,8) has a lower energy than conflguration (2,7) 
for low applied magnetic field. For increasing magnetic 
field, the reverse is true. The rearrangement of the vor- 
tex configurations from the state (1,8) to the state (2,7) 
is related to the change in the steepness of the poten- 
tial energy profile (i.e., the vortex-surface interaction, see 
Eq. ([2)) for different points in the R — H phase diagram. 
The inset of Fig. 7 shows the energy profiles correspond- 
ing to points C, D, and E in the phase diagram. Previ- 
ously, it was shown for charged particles that the particle 
configiration is influenced by the steepness of the conflne- 
ment potential [2^. 



As it follows from the phase diagram (Fig. 7), we can 
expect that for larger radii (in the mesoscopic regime, 
i.e., when R < Xeff', we show in Sec. Ill that for thick 
disks with R > Xeff the conflguration (2,7) restores as 
the ground-state conflguration). the conflguration (1,8) 
has a lower energy than the state (2,7), in the low mag- 
netic fleld range. However, the difference in energy be- 
tween the states (2,7) and (1,8) decreases for increasing 
R, although the energy of the state (1,8) remains lower 
than that of the state (2,7) for even larger disks (e.g., 
with radius i? ~ 40 — 90^ but still "mesoscopic" due to 
the condition R < Xeff valid for very thin disks). This 
seems to contradict our previous result that the state 
(2,7) is the highest probability state in lar ge d isks (which 
is also in agreement with the experiment 'l0|). Thus the 
question is, whether the highest probability conflguration 
(e.g., state (2,7)) is always the ground state? If not, what 
is the reason for that? 

To answer this question, we calculated, using the sta- 
tistical approach described above, the R — H diagram for 
the state (2,7) (Fig. 8), i.e., the region of the existence of 
the state (2,7) and the region where the state (2,7) has 
the highest probability, for radii R — 40^ to 90^. The 
region of the existence of the state (2,7) as the highest 
probability state is very narrow although it is well-deflned 
even for very large radii (e.g., R = 90^, see Fig. 8). How- 
ever, for radii R < 40^ this region becomes narrower and 
unstable, i.e., it can even disappear, giving rise to higher 
probability of appearance of the state (1,8). 

This calculation clearly shows that the highest prob- 
ability state (2,7) is not the ground state in large disks. 
The reason of such behavior is that the energy minimum 
in conflgurational space corresponding to the state (2,7) 
is very wide while the competing state (1,8) possesses, 
although slightly deeper, a much narrower minimum. 
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FIG. 8; The R-H diagram for the vortex state (2,7) cal- 
culated using the statistical approach, for large disks with 
radius R — 40^ to 90^. The area between the red curves and 
dots shows the region where the state (2,7) has the maximum 
probability. 

Thus, statistically the system ends up more often into the 
wide minimum corresponding to the state (2,7). This is 
confirmed by a calculation of the potential energy profile 
which is shown in Fig. 9 as a function of the displace- 
ment (i.e., for different vortex configurations, between 
the initial nonequilibrium configuration (2,7) - through 
the equilibrium state (2,7) - and the final state (1,8)) of 
one of the central vortices of the configuration (2,7) dur- 
ing the continuous transition to the configuration (1,8). 
The position of the other vortices is determined by min- 
imizing the energy. The corresponding changes of the 
vortex configuration are shown in the inset of Fig. 9. We 
started from an out-of-equilibrium (2,7) configuration, 
passed through the equilibrium (2,7) configuration (wide 
minimum), then passed over the energy maximum, and, 
finally, ended up at the equilibrium (1,8) configuration (a 
tail of the transition is also shown for out-of-equilibrium 
(1,8) configuration) which has slightly lower energy. 



III. MACROSCOPIC REGIME: THICK DISKS 

In Sec. II we analyzed in detail the formation of vor- 
tex shells in mesoscopic disks. Although we considered 
rather large disks with radii up to i? ~ 100^, the results 
obtained in the previous section refer essentially to thin 
disks: only if the disk's thickness d is small enough, the 
R < Xeff — X^/d condition is satisfied, i.e., for disks 
with R ~ 1/im, d has to be of the order of a few nanome- 
ters. Nb disks used in the experiment [16] had radius 
i? = 1 — 2/im and thickness d = 150nm. For such thick 
disks with R > X (e.g., A(0) — 90nm for Nb disks in 
the experiment [3 ) : the effects due to London screening 
become important. In this section, we consider the limit 
of thick disks with d 3> A and we study how the Lon- 
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FIG. 9: The energy of the system of L = 9 vortices as a 
function of the displacement rd (measured in units of the disk 
radius rg and counted from the center of the disk) of one 
of the central vortices of the initial configuration (2,7) during 
the continuous transition to the configuration (1,8). The inset 
shows the corresponding evolution of the vortex configuration 
(different symbols show vortex configurations for different dis- 
palcements of the vortex marked by red dashed line) during 
the transition from the state (2,7) to the state (1,8). 

don screening in the vortex-vortex and vortex-boundary 
interactions influences the vortex patterns in the disk. 

Here we model a cylinder with radius R infinitely long 
in the z-direction by a two-dimensional (2D) (in the xy- 
plane) disk, assuming the vortex lines are parallel to the 
cylinder axis. This approach was used for studying, e.g., 
vortex dynamics in periodic [1^ Isol. 3l| and quasiperi- 
odic arrays of pinning sites (APS) 3^. As distinct from 
infinite APSs where periodic boundary conditions are im- 
posed at the boundaries of a simulation cell, here we 
impose boundary conditions at the edge of the (finite- 
size) disk, namely, the potential barrier for vortex en- 
try/exit. To study the configurations of vortices inter- 
acting with each other and with the potential barrier, we 
perform simulated annealing simulations by numerically 
integrating the overdamped equations of motion (unlike 
in Refs. |2^, Jilj, ^52] , there is no external driving force 
in our system, and we study the relaxation of initially 
randomly distributed vortices to the ground-state vortex 
configuration): 



pvb 



(10) 



Here, is the total force per unit length acting on vor- 
tex i, f™ and f"*" are the forces due to vortex- vortex 
and vortex-barrier interactions, respectively, and is 
the thermal stochastic force; rj is the viscosity, which is 
set here to unity. The force due to the interaction of the 
i-th vortex with other vortices is 



j 



(11) 
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where Ny is the number of vortices, Ki is the modified 
Bessel function, f ^ = (r^ — ) / | r j — | , and 

" 87r2A3 ■ 

It is convenient, following the notation used in Refs. 0, 
[boI . [Sll . [3^ , to express now all the lengths in units of A 
and all the fields in units of ^o/A^. The Bessel function 
Ki(r) decays exponentially for r larger than A, there- 
fore it is safe to cut off the (negligible) force for distances 
larger than 5A. In our calculations, the logarithmic diver- 
gence of the vortex- vortex interaction forces for r — > is 
eliminated by using a cutoff for distances less than O.IA. 

Vortex interaction with the edge is modelled by im- 
plying the usual Bean-Levingston barrier 0,[2l, [s^. We 
assume that the repulsive force exerted by the surface 
current on the vortex at a distance r from the disk edge 
decays as 



cvbc 



47rA 



-cxp 



(12) 



as it does in the case of a semi-infinite superconductor 
0) IHj HI] (which is justified for disks with R:s> X), and 
the attractive force due to the vortex interaction with its 
image is expressed by 



and 



cvb 



pvbc 



pvbi 



(13) 



(14) 



Here we assume that for large enough disks, the distance 
from the edge to the image is equal to the distance to the 
vortex. 

The temperature contribution to Eq. (jlOp is repre- 
sented by a stochastic term obeying the following con- 
ditions: 



and 



(/f W) = 0, (15) 



{fnt)fnt'))=2r^ksTS,,S{t-t'). (16) 



The ground state of a system of vortices is obtained as 
follows. First, we set a high value for the temperature, 
to let vortices move randomly. Then, the temperature is 
gradually decreased down to T = 0. When cooling down, 
vortices interacting with each other and with the edges 
adjust themselves to minimize the energy, simulating the 
field-cooled experiments (see, e.g., fs^-lssj). 

Our calculations show that most of the vortex config- 
urations found in Sec. II and in the previous theoretical 
works on mesoscopic disks [26| remain unchanged also 
in large disks where the interactions are screened at the 
London penetration depth A. These are stable shell pat- 
terns (e.g., (1,6) for L = 7 and (1,7) for L = 8, etc.) 
which were found to be the ground-state configurations 
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FIG. 10: The R — H phase diagram for vortex states with 
total vorticity L = 9 in case of a thick disk. The area between 
the dashed blue curves shows the region where the states with 
L = 9 are the ground state. The boundary separating the 
states (1,8) (small radii) and (2,7) (large radii) are shown by 
solid red squares. The insets show the corresponding vortex 
patterns, (1,8) and (2,7). 

of vortices in superconductors [26l] , in liquid He [l3l , and 
in a system of charged particles confined by a parabolic 
potential . These stable configurations are mainly de- 
termined by the circular shape of the disk, and they are 
to a much less extend sensitive to the specific interac- 
tion potentials between the particles and the boundaries. 
On the other hand, the "borderline" configurations (i.e., 
those for which one or more shells start to be filled), e.g., 
the states (1,8) vs. (2,7) for L = 9 or the states (2,9) vs. 
(3,8) for L = II, are much more sensitive to the inter- 
actions in the disk. For example, ior L — 11 the theory 
predicts the configuration (3,8) to be the ground state for 
vortices in He [ij] and for charged particles [l^l, and it 
is was also observed in the experiment [l6j in large disks, 
while the theory predicted [26| the configuration (2,9) in 
small mesoscopic disks. For L — 9 the theory predicts 
that the configuration (1,8) is the ground state for vor- 
tices in He jl7j | and in small mesoscopic superconducting 
disks [2^, while for charged particles the configuration 
(2,7) was predicted (23|. This vortex configuration, (2,7), 
was also observed in the experiment ^16,] with Nb disks. 
In Sec. II we showed that in mesoscopic disks the state 
(2,7) had the highest probability to appear (due to the 
wide potential energy minimum related to this state) , al- 
though it was not the lowest energy configuration, but 
instead the (1,8) configuration was the ground state. 



10 



In case of thick disks as considered here, the calcula- 
tions show the crossover behavior of the vortex patterns 
from the state (1,8) to (2,7) with increasing radius of the 
disk. The phase diagram in Fig. 10 illustrates this be- 
havior. Note that the potential barrier at the disk edge 
becomes extremely low for low values of the applied mag- 
netic field Hq when we have a large radius of the disk. 
This means that it is very difhcult to stabilize a vortex 
state with only few vortices in such a large disk (in ex- 
periment, and also in the numerical calculations using 
the Ginzburg-Landau equations) because for even very 
low barrier at the boundary many vortices can enter the 
sample without any appreciable change of the flux inside 
the disk. The lines separating the states with different 
vorticities (shown by dashed lines in Fig. 10) are calcu- 
lated here assuming the flux inside the disk is on average 
equal to the applied magnetic field Hq multiplied by the 
area of the disk. The calculated line separating the states 
(1,8) and (2,7) is shown by solid squares. The phase di- 
agram shows that in relatively small disks with radius 
i? < 5A (that is R < 30^ in case of Nb disks [ll|), the 
configuration (1,8) is the ground state while for larger 
disks we find the state (2,7), in agreement with the ex- 
periment [l6|. 

This crossover behaviour could be understood in the 
following way. In Fig. 11(a) we plot the vortex confine- 
ment potential profiles for a mesoscopic disk {R ~ A) and 
for a macroscopic disk {R > X). In a mesoscopic disk, 
all the vortices interact with the screening current which 
extends inside the disk. In a macroscopic disk, only the 
outer-shell vortices feel the screening current. More im- 
portantly, the intervortcx interaction changes in a disk 
with the London screening: in a mesoscopic disk, each 
vortex interacts with all other vortices since the currents 
created by the vortices strongly overlap (see Fig. 11(b)), 
and the minimum potential energy is reached when the 
sum of all the intervortex distances is maximum, i.e., for 
the configuration (1,8). In a macroscopic disk, the inter- 
vortex interaction is very weak, and each vortex interacts 
only with its closest neighbor through the tails of the cur- 
rents associated with each vortex (see Fig. 11(c)), and the 
minimum potential energy is reached when the sum of 
closest-neighbor intervortex distances is maximum, i.e., 
for the configuration (2,7). The vortex pattern (2,7) in 
a large disk (see inset in Fig. 10) resembles a distorted 
Abrikosov vortex lattice in an infinite superconductor 
which is stabilized by intervortex interactions in the ab- 
sence of boundaries (note that the outer-shell vortices are 
relatively closer to the boudary and the two vortices in 
the inner shell are slightly out of the center minimizing 
the interaction energy with the 2 and 3 neighbors). 

The calculated crossover behavior found here is consis- 
tent with the R — H phase diagram obtained in Sec. II 
for mesoscopic disks (Fig. 7) that predicted the config- 
uration of (1,8) as the ground state for radii R > 10^. 
Thus, according to the phase diagrams for mesoscopic 
disks (Fig. 7) and for macroscopic disks (Fig. 10), there 
are two crossovers between the states (1,8) and (2,7): the 
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FIG. 11: The vortex confinement energy profiles for: meso- 
scopic disk {R ~ A) and for macroscopic disk [R > A) (a). 
Schematic plots illustrating the magnetic field profiles (green 
dashed curves) of the interacting vortices in: mesoscopic disk 
(solid red curve shows the vortex confinement potential) (b), 
and in macroscopic disk (solid blue curve shows the vortex 
confinement potential) (c). 

configuration (1,8) is the ground state in disks with ra- 
dius 10^ ^ R ^ 30^, while the configuration (2,7) occurs 
to be the ground state in large disks with R > 5X (i.e., 
R > 30^ in Nb disks [lB|) and in very small disks with 
R < 10^. The mechanism of the second crossover for very 
small disks is very different from that for large disks, and 
the transition (l,8)-^(2,7) happens in very small disks 
due to a strong overlap of the vortex cores in the outer 
eight-vortex shell: the vortices cannot accomodate on the 
outer shell and one of them is pushed towards the interior 
of the disk (note that for even smaller disks the configura- 
tion (2,7) collapses to a giant-vortex state). This behav- 
ior will be demonstrated in Sec. IV using the Ginzburg- 
Landau theory. 



IV. COMPARISON WITH THE GL THEORY 

In order to go beyond the London approximation we 
used also the Ginzburg-Landau (GL) equations to calcu- 
late the free energy and find the ground state. Within the 
GL approach vortices are no longer point-like "particles" 
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FIG. 12: Contour plots of the Cooper-pair density of state 

(2.7) for R = 8.3^ (a), and state (1,8) for 7? = 8.2^ (b), in an 
applied magnetic field Ho = 0.35-ffc2- Blue and red regions 
correspond to low and high Cooper-pair density, respectively. 

but extended objects. The expression for the dimension- 
less Gibbs free energy is (see, e.g., Ref. [4]): 

G = V-^ [ [2(A - Ao) • ho- I * |^]dr, (17) 
Jv 

with iI}{t) the order parameter, A(Ao) the vector poten- 
tial of the total (applied) magnetic field and j2D the su- 
perconducting current. By comparing the dimensionless 
Gibbs free energies of the different vortex configurations, 
we find the ground state. Similarly we could find the 
two stable configurations (2,7) and (1,8) in a disk with 
vorticity L = 9 as we found within the MD simulations 
in Sees II and III. 

The results of our calculations of the order parameter 
distribution using the GL equations (for simplicity, this 
was done for zero disk thickness, c? — > 0, i.e., in the limit 
of extreme type II superconductor, for a given applied 
magnetic field, i.e., only first GL equation was solved) 
for the total vorticity L = 9 are shown in Fig. 12. The 
states (2,7) and (1,8) are shown in the phase diagram (see 
Fig. 7) by symbols A and B, respectively. Samples witli 
different radius were examined for a fixed external mag- 
netic field H = 0.35Hc2- For a disk with radius R — 8.2^, 
our calculation gives (1,8) as the ground vortex state. 
When the radius of the disk is increased, the energy of 

(1.8) is the lowest one till R — 8.25^, after which config- 
uration (2,7) becomes the ground state. This compares 
with our results of Sec. II using the London theory where 
we found that the transition (1,8) (2,7) occurred for 
R = 8.13^ when H = 0.35i?c2- Thus, the results of a cal- 
culation of the vortex configurations within the GL, with 
an appropriate choice of the radius and external param- 
eters, confirms the crossover behavior found in Sec. II. 

V. CONCLUSIONS 

In this work, we studied the vortex configurations 
in mesoscopic superconducting thin disks and in thick 
disks taking into account the London screening, using 
the Molecular-Dynamics simulations of the Langevin- 
type equations of motion and confirmed these results, in 



case of small disks, using the more extended Ginzburg- 
Landau functional theory. 

This study was motivated by recent experiments by 
Grigorieva et al. who observed vortex shell struc- 

tures in mesoscopic Nb disks with i? ~ 1 — 2.5/im by 
means of the Bitter decoration technique. It was shown 
in those experiments, that in disks with vorticity rang- 
ing from L = 1 to 40, vortices fill the disk according 
to specific rules, forming well-defined shell structures, as 
earlier predicted in Ref. J22]. They analyzed the forma- 
tion of these shells which resulted in a "periodic table" of 
formation of shells. It was shown that most of the experi- 
mentally observed configurations for small L agreed with 
those theoretically predicted earlier [l^, |2^ . At the same 
time, some of the configurations which were observed in 
these experiments were not found earlier in vortex sys- 
tems (although they were shown to appear in systems of 
charged particles and in superfluids). 

In this work, we found the rules according to which the 
shells are filled with vortices for increasing applied mag- 
netic field. In particular, it was shown in our calculations, 
that for the vortex configurations with the number of vor- 
tices up to L = 5, the vortices form a single shell. The 
formation of a second shell starts from L = 6. Similarly, 
the formation of a third shell starts at L = 17, and of a 
fourth shell at L = 33. These theoretical findings are in 
agreement with the results of the experimental observa- 
tions of Ref. [l6|. Moreover, we found those states which 
were observed in the experiments but not found in pre- 
vious calculations. Thus, we filled the missing states in 
the "periodic table" of vortex shells in mesoscopic disks. 
We studied in detail the region of parameter space where 
those states exist, and compared the obtained results to 
previous theoretical works where small mesoscopic disks 
with iZ ~ 5 — 10^ were considered. 

It was shown that some of the vortex configurations 
(i.e., those which are at the borderline between configu- 
rations characterized by different stable shell structures) 
are very sensitive to the size of the disk. For instance, 
we found that depending on the radius of the disk, there 
are two crossovers between the states (1,8) and (2,9) for 
L = 9: at i? ~ 10^ and R ~ 30^. The (1,8) (2,7) tran- 
sition occurs for disks with R ^ 5X (that corresponds to 
R ^ 30^ in case of Nb disks in the experiment [16] ) due 
to the effect of the London screening in large disks, while 
in small disks with R ^ 10^ this transition happens due 
to the compression of the outer eight-vortex shell. 

Thus we performed a systematic study of the size- 
dependence of vortex configurations in mesoscopic su- 
perconducting disks. Our results agree with the experi- 
mental observations of vortex shells in Nb disks [iGj and 
explain the revealed discrepancies with the earlier calcu- 
lations of vortex shells. 
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